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A  Fortran  Code  for  Calculation  of  Eigenvalues  and 
Eigenfunctions  in  Real  Potential  Wells 

Randall  S.  Caswell 


A  Fortran  code  has  been  developed  for  the  calcu- 
lation of  eigenvalues  and  eigenfunctions  for  neutrons 
in  real  potential  wells.  A  systematic  procedure  is 
given  for  approximate  location  of  the  eigenvalue  and 
an  automatic  search  procedure  to  determine  the  exact 
location.   The  code  may  be  used  for  either  bound  or 
scattering  states.   In  the  case  of  scattering  states, 
the  criterion  for  maximum  scattering  (90°  phase  shift) 
is  used  to  determine  the  energy  of  the  state.   The 
eigenvalues  are  determined  by  matching  the  numerically 
calculated  logarithmic  derivative  (f  .)  inside  the  nucleus 
to  the  appropriate  analytical  logarithmic  derivative  for 
the  region  outside  the  nucleus.   In  an  alternate  mode  of 
operation,  the  outside  value  of  f .  may  be  set  arbitrarily, 
and  a  match  made  to  this  value.  Sample  results  for  a 
Woods-Saxon  well  with  spin-orbit  coupling  for  the  case 
of  oxygen-16  are  shown.  The  code  is  in  Fortran  and  was 
written  for  an  IBM  7090  computer. 

1.   Introduction 

Nuclear  shell  model  calculations  have  until  now  employed  almost 
exclusively  harmonic  oscillator  wave  functions  for  the  evaluation  of 
radial  integrals.  This  has  been  done  both  for  bound  states  and  for  states 
which  actually  belong  to  the  continuum.   The  bound  states  evidently  can  be 
very  well  represented  by  harmonic  oscillator  wave  functions.  Any  errors 
due  to  the  difference  between  the  "actual"  optical  potential  and  the  har- 
monic oscillator  potential  will  in  general  be  much  smaller  than  the  errors 
due  to  neglect  of  higher  configurations.   This  is,  however,  not  necessarily 
true  for  continuum  states. 

The  bound  state  wave  functions  have  quite  a  different  character  than 
the  continuum  wave  functions.   In  a  treatment  of  continuum  states  it  is 
thus  desirable  to  have  something  better  to  work  with  than  the  harmonic 

oscillator  eigenfunctions.   The  present  report  is  concerned  with  the 
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evaluation  of  eigenf unctions  specified  by  a  given  potential  well  and 
boundary  conditions. 

Calculations  of  eigenvalues  in  a  diffuse  potential  have  been  re- 
ported by  Ross,  Mark,  and  Lawson  (195&)  where  they  used  an  analog  com- 
puter and  obtained  eigenvalues  accurate  to  about  0.1  Mev.  The  present 
calculation  is  for  a  digital  computer  (IBM  7O9O) ,  is  inherently  capable 
of  much  higher  accuracy,  and  also  provides  numerically  integrated  eigen- 
f unctions  as  required. 

2.  Description  of  the  Calculation 

This  code  has  been  developed  from  a  nuclear  optical  model  code 
written  for  calculation  of  neutron  elastic  scattering  which  has  been 
described  elsewhere  (Caswell,  I962) .  As  a  result  some  of  the  features 
are  slightly  more  complicated  than  would  be  necessary  for  the  present 
purpose o  The  main  calculational  steps  are  as  follows:   (l)  input  of  data; 
(2)  calculation  of  the  real  potential  well,  which  is  a  Saxon  well  with 
spin-orbit  coupling  in  the  present  code,  but  can  easily  be  changed;  (5) 
numerical  integration  of  the  wave  functions,  u(r)  =  rt(r) ,  out  to  a  maxi- 
mum radius  which  is  taken  where  the  nuclear  potential  has  become  negligibly 
small  (except  during  the  search  procedure  where  smaller  radii  are  used), 
and  calculation  of  the  logarithmic  derivative,  f .(inside) ;  (k)    calculation 
of  the  required  f . (outside)  if  it  is  not  fixed  in  the  input  data.  The 
approximate  location  of  the  eigenvalues  is  determined  by  a  survey  versus 

energy  of  f. (inside),  f. (outside),  and  the  difference  between  them.   An 

> 
automatic  search  procedure  finds  the  energy  for  which  f .(inside)  = 

f. (outside),  and  the  wave  function  vs.  radius  is  calculated  and  printed 

or  punched  on  cards  for  the  energy  of  the  eigenvalue. 


The  potential  is  calculated  by  subroutine  VR3  as  a  function  of 
radius  at  up  to  500  points,   Spacing  between  points  is  determined  by  the 
variable  HI  up  to  radius  RI,  by  HII  up  to  RII,  and  by  HIII  up  to  RMAX. 
The  potential  used  is; 

V(r)  -  Vcp(r)  -  aVc  Q^   (1-3)1  M  (1) 

where  p(r)  =  - ,  the  "Saxon"  potential;  a,   is  the  strength 

1  +  ***<rTV 

of  the  spin-orbit  interaction  compared  to  the  "Thomas"  interaction  for  a 
nucleon;  M  is  the  nucleon  mass;  V  is  the  central  real  potential;  h/Mc  is 

the  Compton  wave  length  for  a  nucleon;  i   is  the  orbital  angular  momentum 

— * 
of  the  neutron;  c   is  the  Pauli  spin  operator  of  the  incident  neutron;  R 

is  the  nuclear  radius;  and  a  is  the  "dif fuseness"  parameter  of  the  poten- 
tial.  If  j  =  I   +  Jj  (chosen  by  letting  variable  MSPIN  =  l)  ,  then  o'i  =  I, 
and  the  values  of  the  variable  VPLUS  (potential  for  spin-orbit  parallel) 
for  different  radii  are  calculated.   If  j  =  i   -  ■§•,   (chosen  by  letting 
MSPIN  =  -1),  then  o-i  =   -(i+l) ,  and  the  VMINUS  values  (potential  for  spin- 
orbit  anti-parallel)  are  calculated.   Other  potential  shapes  could  be  used 
by  simple  changes  in  this  subroutine.   The  four  arbitrary  parameters,  VC, 
ALPHA,  RADIUS,  and  A  could  be  used  to  define  the  new  potential  with  no 
changes  outside  of  subroutine  VR3.  Any  arbitrary  potential  could  be  read 
in  point  by  point  by  replacing  VR3  with  an  appropriate  "READ"  statement. 

The  calculation  of  f. (inside)  and  f. (outside),  and  the  search  proce- 
dure for  matching  the  values  is  controlled  by  subroutine  MATCH.  First 
the  radial  part  of  the  non-relativistic  Schr$dinger  wave  equation  is 


integrated  out  in  radius  for  the  appropriate  i\ 

o 
d  U/r)   2m  f_   lU+l)h2  „,   vl   ,  v    _         ,. 

T^~~  ^  LE  "  t  2   "  V(r) Jui(r)  =  °"       (2) 

dr       h        2mr 
The  values  of  the  function  and  the  first  derivative  at  the  first  point 
(radius  =  Hi)  are  determined  from  the  relation  u.(r)  =  C.r    as  r  -*  0, 
where  C.  is  a  constant.   The  second  order  Runge-Kutta  method  is  used  for 
the  integration  (see  Zurmtihl  196l)  .   The  integration  is  carried  out  by- 
subroutine  INTEG3  and  is  stopped  when  the  radius  exceeds  RSTOP,  and 
f.  (inside)  =  R(du./dr)/u.  is  calculated., 

Calculation  of  f  (outside) .   The  calculation  may  be  run  in  two  alter- 
native modes.  When  MODE  =  1,  if  the  energy  is  negative  (bound  state),  the 
outside  wave  function  is  required  to  have  an  f .  which  corresponds  to  an 
exponential  function  which  asymptotically  approaches  zero  for  large  r. 
If  the  energy  is  positive  (scattering  state),  the  criterion  of  a  90° 
phase  shift  is  applied  to  the  wave  function  at  the  radius  RSTOP  where  the 
inside  and  outside  functions  must  join  smoothly.  When  MODE  =  2,  then  f . 
is  read  in  from  the  input  data  cards  and  no  calculation  of  f.   is  made. 

The  outside  wave  function  is  the  solution  of  equation  (2)  in  the 
region  of  space  where  the  nuclear  potential  V(r)  =  0,  but  the  centrifugal 
potential  term  is  present.   The  solutions  in  general  are  spherical  Bessel 
functions.   For  bound  states,  the  "outside"  or  asymptotic  solutions  are 
spherical  Hankel  functions, 


^j&M  =  LXVX)  +  ixni(x) j  =  F£  +  iGi  where  x  =  ^Ajt  lEl  r* 

For  the  90°  phase  shift  continuum  case,  the  solutions  are  simply 

xn^(x)  =  G^.   However  for  matching  we  require  not  the  function  but  the 
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logarithmic  derivative,  f .  <,      It  is  possible  to  calculate  f. (outside)  in  a 
number  of  equivalent  ways,  for  example  by  first  evaluating  the  functions 
and  then  the  derivatives  from  appropriate  recursion  relations.   The  cal- 
culation as  done  here  uses  the  quotient  functions  of  the  Bessel  functions 
discussed  by  Onoe  (1958)  f°r  bound  states,  and  a  relation  given  by  Blatt 
and  Weisskopf  (1952)  for  the  continuum  case.  Details  are  given  in  Appen- 
dix 1. 

3.  Search  Procedure 
If  the  search  for  eigenvalues  is  carried  out  by  matching  the  f  *s  at 
some  distance  outside  the  nuclear  potential,  particularly  for  states 
located  deep  in  the  potential  well,  there  is  a  high  probability  of  miss- 
ing the  state  completely.  The  reason  for  this  is  that  during  the  numerical 
integration  of  the  wave  function,  as  soon  as  the  first  radial  point  out- 
side the  nuclear  potential  at  that  energy  is  reached,  the  eigenf unction 
becomes  in  general  a  rapidly  growing  exponential,  and  f .  becomes  very 
large.  Only  when  the  energy  is  chosen  extremely  close  to  the  correct 
eigenvalue  does  the  wave  function  taper  toward  zero  and  have  a  small  or 
negative  value  of  f..  The  quantities  f.  are  thus  very  sensitive  functions 
of  the  energy.   In  order  to  systematically  find  every  eigenvalue,  some  pro- 
cedure is  needed  to  survey  the  eigenvalues  under  conditions  of  low  sensi- 
tivity. The  calculation  was  therefore  split  into  two  parts:   first  we 
perform  a  survey,  and  then  we  use  an  automatic  search  procedure  starting 
with  relatively  small  matching  radius  (which  means  low  sensitivity)  and 
increasing  it  stepwise  to  the  desired  value.  This  requires  two  separate 
runs  on  the  computer,  the  first  for  the  survey  procedure  to  find  the 
approximate  eigenvalue  location,  and  the  second  to  find  the  exact  location 
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by  the  automatic  search  procedure. 

Survey  procedure.  When  the  variable  NOMAX  =  1,  only  one  integration 
out  in  radius  is  made  at  each  energy,  and  that  integration  is  stopped  at 
the  edge  of  the  well  (for  the  particular  energy)  for  negative  energy,  and 
at  the  value  of  RSTART  for  positive  energy  or  for  negative  energy  if  RSTART 
corresponds  to  a  smaller  radius  than  the  first  criterion  above.  The  maxi- 
mum radii  of  integration  for  a  typical  situation  are  shown  as  the  heavy 
black  line  in  figure  1.  The  values  of  f. (inside),  f. (outside)  and  the 
difference  between  them  are  printed  out  at  the  specified  input  energies 
taken  typically  0.5  or  1  Mev  apart.  Typical  results  of  this  survey, 
called  the  "signature  of  the  well"  are  shown  in  figures  2  and  3  f°r  the 
case  of  si  and  p~/_  potential  wells  respectively.  The  small  arrows 
indicate  the  accurate  eigenvalue  locations  found  by  the  automatic  search 
procedure.  The  approximate  eigenvalue  locations  corresponding  to  differ- 
ence =  0,  are  used  as  starting  energies  for  the  automatic  search  procedure. 

Automatic  search  procedure.  First  an  integration  of  the  wave  func- 
tion is  carried  out  at  the  approximate  eigenvalue  energy.  A  second  inte- 
gration is  carried  out  at  a  nearby  energy.  Using  the  differences  (between 
f  *s)  found  in  these  two  trial  calculations,  a  linear  interpolation  (or 
extrapolation)  is  made  to  an  energy  for  which  the  difference  is  predicted 
to  be  zero.  This  is  the  predicted  eigenvalue  energy.   Integration  is 
carried  out  at  the  third  energy,  and  the  difference  is  again  found.  We 
now  know  three  energies  and  the  three  corresponding  differences .   Based 
on  this  information,  by  quadratic  interpolation,  we  predict  the  fourth 
energy  for  which  the  difference  should  be  closer  to  zero.  The  calculation 
continues,  energy  prediction  being  made  by  quadratic  interpolation,  until 
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either  the  difference  is  closer  to  zero  than  the  difference  specified  in 
the  input  by  the  variable  EPS,  or  until  the  maximum  number  of  tries  has 
been  used  up.   The  maximum  number  of  tries  is  determined  by  the  variable 
NOMAX. 

The  computer  time  for  the  sample  automatic  search  procedure  shown  in 
Appendix  kt   which  required  21  repetitions  of  the  integration  of  the  wave 
function  and  matching  procedure  was  0,5  min.   For  a  systematic  search  pro- 
cedure at  69  energies,  the  computer  time  was  O.k   min. 

The  matching  described  above  is  carried  out  at  a  radius  determined  by 
the  variable  RSTART.   RSTART  should  be  chosen  in  the  region  where  sensi- 
tivity of  the  f .   changes  with  energy  is  low.  Having  determined  the  eigen- 
value  for  a  potential  which  is  cut  off  at  a  value  of  radius  determined  by 
RSTART,  we  use  this  energy  as  the  first  trial  energy  for  a  larger  value  of 
r,  which  is  determined  by  the  variable  RSTOP.  The  permissible  size  of  a 
step  in  radius  may  be  estimated  from  the  WKB  approximation.   If  we  let 
cp  =  f  (inside)  -  f. (outside),  then  let  S(R)  =  -r^,  a  quantity  which  we 
have  previously  called  the  sensitivity.   In  the  outward  steps  in  radius, 


<Rn+l) 


the  increase  in  sensitivity,  — >'  ^  =  c  ,  should  be  kept  about  constant 

^  n' 

for  each  increase  in  RSTOP.  Approximately  we  may  obtain  o_  as  follows: 

K 

R     .    .      ^_^_ 

£n  ctr  =    f  £  72m(V-E)    dr 

R 
n 


R 

a4 


Jn+1£y^[I!    dr^/^AR 


R 
n 

since  V  soon  becomes  small  in  this  region  outside  the  well.  The  change  in 

4nCJ 
RSTOP  is  given  by 


AR~y2mM7£ 

10 


where  0"  is  an  arbitrary  parameter  determine i  by  the  input. 
R 

At  the  new  matching  radius,  the  first   stimate  of  energy  is  the 

converged  eigenvalue  at  the  previous  matching  radius.  The  second  estimate 

2 
should  be  lower  since  VR  for  the  well  will  be  slightly  larger,  causing 

the  new  eigenvalue  to  be  slightly  lower.  At  the  second  matching  radius, 

the  second  trial  eigenvalue  is  generated  in  the  same  way  as  at  the  first 

matching  radius.  At  third  and  higher  matching  radii  it  is  generated  from 

the  relation: 

AE  .  =  i-  (AE  ) 
n+1   a„  rr 

£1 

where  o_  is  an  arbitrary  parameter  determined  from  the  input.  Figure  k 
shows  a  typical  case  of  converging  on  an  eigenvalue  using  this  procedure 
of  calculating  at  successive  matching  radii.   The  final  match  is  carried 
out  at  the  radius  RMAX.  When  a  converged  eigenvalue  is  found  at  radius 

RMAX,  the  integration  in  radius  is  repeated  and  the  wave  function  vs. 

f*  2 
radius  is  printed  out.  Also  the  val    of  thi    u  dr  is  calculated  by  the 

trapezoidal  rule  for  normalization  purposes,  and  is  printed  out. 

ko      Sample  Results 

Sample  results  of  the  automatic  search  procedure  for  the  case  of 

0xygen-l6  are  shown  in  figure  h,      -' "'-?envalues  are  shown  for  the  Is,  2s, 

3s,  the  two  lp  states,  the  Id  state,  and  the  eigenf unctions  are  shown  for 

the  three  s  states. 
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Appendix  1.   Evaluation  of  f  (outside) . 
For  bound  states  the  logarithmic  derivative,  f  ,  may  be  obtained 
from  one  of  the  quotient  functions  of  the  Bessel  functions  which  have 


been  discussed  by  Onoe  (1958) . 

,-du 

f, (outside)    =  R-^^ 


Cdr  )  R  __       Ldr   (rhi). 


=   X 


(UJL)R  [r  hjt 


R 


x  (3) 


hx(x)jx 


where  X  =  i  /^|e|r. 

but  dT  <x  V  =  k  G*+1*"  Vx0  -  (*+<Vx0  -  x  Wx)  •       w 


■'  '  A 


Therefore. 

'l+l   -  X  h     X(X)  )                    X  hjtfl(x) 
^(outside)    =  X  T^-(T) =  Ul   -       ^x) (5) 

For  a  continuum  state  we  have  a  similar  relation  in  terms  of  the 
Neumann  functions  rather  than  the  Hankel  functions: 

x  n»+i(x)  /-     G.\ 

h  -  »\- -$r(r  * £>  (6) 

where  the  prime  refers  to  differentiation  with  respect  to  X.  For  a 

bound  state,  /,\ 

^(outside)  =  i+1  -^t(x), 

in  the  notation  of  Onoe. 

In  the  notation  of  Magnus  and  Oberhettinger  (195M» 

x  Ku+1(x) 
C\     (l)  =   >  \ — — .  A  recursion  relation  may  be  derived  which 


Ik 


~  (1) 

is  used  by  the  computer  to  calculate  x/v  '  for  the  appropriate  value  of 


I: 


^U+  x 


w. 


frv-l 


(7) 


For  continuum  states,  the  criterion  of  a  90*  phase  shift,  corre- 
sponding to  a  maximum  value  of  the  scattering  cross  section  (or  "resonance") 
is  used  to  determine  the  location  of  the  state.  The  relation  for  f .   may 
be  written  as  in  equation  (6)  above 


or 


h  =  Az  +  h  tan  h> 


(8) 


V  i  ¥i  rG4G/+F// 

where  exp(-2i^)    =  ^   .  ±  ^  =  \G ^  +  F^2     _ 


r=R; 


and  S„   =  R 


Vi'  -  fa' 


1    XV  +  V 


r=R.   (See  Blatt  and  Weisskopf ,  1952) 


Substituting  for  exp(-2i§.)  we  have 


h  -  h  -  h 


<*tGt 


/Gi2 "  F/> +  <F/ +  G/>] 


(9) 


Equation  (9)  is  equivalent  to  equation  (5)»  and  was  the  formula  actually 
used  in  the  calculation.  Analytic  expressions  for  all  of  the  quantities 
in  equation  (9)  appear  in  the  report  of  Amster  and  Culpepper  (1957)* 
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Appendix  2.  Meaning  of  the  Input  Variables 

The  meaning  of  the  input  variables  are  listed  here  for  convenience 

in  using  the  code.   Format  may  be  obtained  from  Appendix  5* 

RADIUS     radius  in  the  Saxon  potential  in  fermis 

A         diffuseness  parameter  of  the  Saxon  potential  in  fermis 

VC        central  real  potential  in  Mev  (is  negative  for  attractive 
potentials) 

ALPHA      magnitude  of  the  spin-orbit  potential  expressed  in  number  of 
times  larger  than  for  the  Thomas  term  of  a  nucleon 

AMASS      atomic  mass  in  a.m.u.  of  the  nucleus  except  the  single  neutron 
whose  states  are  being  considered 

CONV       parameter  which  determines  energy  spacing  between  first  two 
trial  eigenvalues,  typical  value  0.010 

MDELTA     the  radial  potential  is  printed  out  at  every  MDELTA'th  point 

JPRINT     conditional  print  variable,  =  1  for  normal  operation,  =  3 
maximum  print  out  for  code  checking 

L         orbital  angular  momentum 

EPS        if  the  difference  between  f  's  is  less  than  EPS,  then  the 
eigenvalue  is  considered  converged. 

MSPIN      =  1  for  spin-orbit  parallel,  =  -1  for  spin-orbit  anti-parallel 

HI        spacing  between  radial  integration  points  at  radii  less  than 
RI,  in  fermis 

HII        same,  but  for  radii  between  RI  and  RII 

HIII       same,  but  for  radii  between  RII  and  RMAX 

RMAX       maximum  radius,  maximum  value  of  RSTOP 

NOMAX      maximum  number  of  tries  for  convergence  on  an  eigenvalue  at  a 

given  value  of  RSTOP,  =  1  for  survey  procedure)  maximum  value  20 

JPUNCH     =  1  for  no  punch  out  of  wave  function,  =  2  for  wave  function 
to  be  punched  on  cards 

KZ        number  of  trial  energies  for  the  given  run,  maximum  =  100. 
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RSTART     initial  value  of  RSTOP 

SIGMAR     variable  which  indirectly  determines  spacing  between  successive 
values  of  radius  (RSTOP)  at  which  eigenvalue  search  is  carried 
out,  typically  2.0 

SIGMAE     variable  which  determines  trial  eigenvalue  spacings  after  the 
first  two  trial  radii,  typically  10.0 

MODE       =  1  for  normal  search  for  bound  and  scattering  states 

=  2  for  use  with  arbitrary  f .   determined  from  the  input 


trial  energy  in  Mev  (negative  for  bound  states,  positive  for 
scattering  states) 


EN(KY) 

FLOuT(KY)   value  of  arbitrary  f  (outside) 
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Appendix  3*   Format  of  Data  Cards 

First  card;    C  in  column  1,  any  desired  information  describing  the  run 
(such  as  name,  date,  element,  etc.)  in  columns  7-72. 


Second  card: 


Third  card: 


Variable 

End 

s  in  Column 

Example 

RADIUS 

8 

3.150 

A 

16 

O.65O 

VC 

2k 

-50.000 

ALPHA 

32 

35.000 

AMASS 

ko 

15.000 

CONV 

kQ 

0.010 

MDELTA 

52 

10 

JPRINT 

56 

1 

L 

60 

0 

EPS 

70      0 

.0100000 

MSPIN 

72 

1 

Variable 

End 

s  in  Column 

Example 

HI 

8 

0.020 

RI 

16 

0.100 

HII 

2k 

0.050 

RII 

32 

1.000 

HIII 

ko 

0.010 

RMAX 

kQ 

6.4oo 

NOMAX 

52 

10 

JPUNCH 

56 

2 

(Note:   decimal  point 
location  may  be  changed 
from  that  shown  here, 
however  print  in  the 
output  will  be  as  shown) 


10   (Note  that  this  and  other 
variables  without  a 
decimal  point  are  inte- 
gers, decimal  point  may 
not  be  used) . 
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Fourth  card:   Variable   Ends  in  Column   Example 


KZ 

4 

1 

RSTART 

12 

4.0 

SIGMAR 

20 

2.0 

SIGMAE 

28 

10.0 

MODE 

32 

1 

Fifth  and  higher  cards:   If  MODE  =  1,  list  the  energies  with  the  last 
character  successively  in  column  8,  16,  24,  32,  40,  48, 
56,  64,  725  using  as  many  cards  as  required  for  the  KZ 
number  of  energies . 

ExamPle:      -29.000  -29.500 

If  MODE  =  2,  list  the  energies,  one  per  card  as  follows: 

Variable   Ends  in  Column   Example 

EN  8  -29,000 

FLOUT  24    0.2000000E+01 

The  variable  FLOUT  is  in  exponential  notsiior  and  in  the  example  means: 

0.2000000  x  10+1  =  2, 
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Appendix  k.     Fortran  Statements  and  Sample  Output 


EIGENVALUES, EIGENFUNCTIONS  JULY    23,1962       0t>230 

READ     10 
10    FORMAT172H 

PRINT     10  J 

DIMENSION    VPLUS(500 )  ,VMINUS(500)  , ENERGY) 10, 100), XP( 1,100), 
1XMP( 1  ,100)  ,XPP  (  1  ,100  )  ,AJLP {  1  -ICO)  ,AJLM( 1 ,100)  ,AILM( 1  ,100)  , 
2AILPI  1  ,100)  ,FL  (  1  ,10  0)  ,XX( 10  0)  ,EN( ICO)  ,SPHTAN( 10,100) ,XM( 1 , 1 00 ) , 
3    REFLt  1,100)  ,DIF(  10     )  ,HT  I  LDE  (  1  » 100  )  ,E  (  1  00  )  ,  PSKA), 

AKRR(A)  ,S( 1  ,100  )  , DELTA ( 1,100),TFG(1»100),FPG(1»100)»GMF(1,100)» 
5P( 10,100) ,G( 10,100) , FLOUT! 100) 

COMMON    VPLUS»VMINUS,XP,XM,XPP,XMP,A-JLP,AJLM,AILP,AILM,FL,RR»X,V, 
1EN, FF, ELL, FRR, RADIUS, HI , A , E , VC » ALPHA ,6 , AMASS ,EMAX ,M ,MDELTA ,R , H I  I , 
2RI  I  ,HI  I  I  ,RMAX,C,N,REXP,VCC,VS,AK,K.Z»KX,KP,KY,K.,XPLUS,XPLUSP, 
3XMINUS,XMINUP,AI  ,AI  I  ,AI I  I  ,CI  ,CI I ,CI  I  I  ,H,XX,F,  I  ,RI  ,MPR I  NT , JPR INT , 
4REFL,SPHTAN,KK,MM,C0NV,N0,EPS,ENERGY-,MSPIN,DIF,HTILDE,N0MAX, 
5JPUNCH,R5T0P,DELTA     ,S IGMAE , S I GMAR ,NTRY , D I FCNV , S , P , PS  I , RRR , S ,GMF ,Q , 
6TFG»FPG»PMOM,MMM,RSTART , FLOUT ,MODE,PROINT ,L 

CALL     INPUT3 

CALL    VR3 

CALL    MATCH 

CALL    ENDJOB 
END 


Note:      Identical  dimension  and  common  statements  must 

be  included  in  every  subroutine,   but  are  omitted 
here  for  brevity. 
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Appendix  4.     Fortran  Statements  and  Sample  Output — continued 

SUBROUTINE     INPUT3 

READ     1 ,RAD1US»A,VC»ALPHA,AMASS»C0NV,MDELTA, JPRINT.L     ,  EPS.  MSP  IN 

PRINT    4501 

4501  FORMATI119H   RADIUS        A       VC   ALPHA     AMASS     CONV  MDL  JP 
1R    L     EPS    MSPIN  ) 

PR  I  NT  1 ,RADIUS,A,VC,ALPHA,AMASS,CONV,MDELTA,JPRINT,L  ,EPS»MSPIN 
1  FORMAT  (6F8.3.3I4,1F10.7»1I2) 
KK=L+1 

READ  7,HI,RI,HII,RII,HIII , RMAX ,NOMAX , JPUNCH 
PRINT  4502 

4502  FORMAT (119H      HI        RI      HI  I      RII     HI  I  I     RMAX  NOMAX 
1JPUNCH  ). 

PRINT  7,HI,RI,HII,RII,HIII , RMAX . NOMAX .JPUNCH 
7  FORMAT(6F8.3.2 14) 
PRINT  4503 

4503  FORMAT(50H  ENGS  RSTART   SIGMAR   SIGMAE   MODE  ) 
READ   4040, KZ, RSTART, SIGMAR, SIGMAE, MODE 

PRINT  4040,  KZ,  RSTART.,  SIGMAR,  SIGMAE,  MODE 
4040  FORMAT! 14, 3F8. 1,14) 

IF  (MODE-1)  4081, 4061, 4082 

4081  PRINT  4504 

4504  FORMAT (50H  TRIAL  ENERGIES  ) 
READ  4080, (EN(KY)  ,KY  =  1,KZ) 

PRINT  4080, (EN(KY) »KY=1 ,KZ  ) 
4080  FORMAT (9F8. 3) 
GO  TO  4085 

4082  PRINT  4083 

4083  FORMAT(50H   ENERGY        FL  OUT  ) 
READ   4084,  ( EN ( KY ), FLOUT ( KY ),  KY  =  1  ,KZ  ) 

PRINT  4084,  ( EM(KY)  ,FLOUT(KY  )  ,KY  =  1  ,KZ  ) 

4084  FORMAT  (F8.3,E16.7) 

4085  CONTINUE 
RETURN 
END 
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Appendix  k.     Fortran  Statements   and  Sample  Output — continued 


SUBROUTINE    VR3 

R  =  HI 

H  =  HI- 

K  =  KK 

M  =  l 

PRINT  2091 
2091  FORMAT! 50H      VPLUS  VMINUS 

2090  DO  2210  1=1 ,1000 

1  =  1 
9  C=  .0110270/A 
15  REXP=EXPF( (R-RADIUS) /A) 

2151  VCC=VC/( 1.0+REXP) 

2152  VS  =  AI_PHA*VC*C*REXP/  (  (  ( 1  .0+REXP  )**2  )  *R-J 
19  AK  =  ,< 

IF  (MSPIN)  22*22.21 

21  VPLUSf  I  )=VCC+VS*(AK-1.0) 
GO  TO  4020 

22  VMINUS ( I )=VCC-VS*AK 
VMINUS! 1 )=VCC 

4020  IF  (M-MDELTA)  3150,28,28 

28  PRINT  30,     VPLUSI  I  )  ,VMINUS(  I  )  ,R 
30  FORMAT(2E15.7,F7.3) 
M  =  0 
3150  M=M+1 
2126  R=R+H 

2130  IF(R-RMAX+.C01 )  2140,2220,2220 
2140  IF(R-RII+.001)  2150,2160»2160 
2150  IF (R-RI+.001 )  2200,2180,2180 
2160  H  =  HI I  I 

GO  TO  2210 
2180  H=HI I 

GO  TO  2210 
2200  H=HI 
2210  CONTINUE 
2220  RETURN 
END 
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Appendix  4„     Fortran  Statements   and  Sample  Output — continued 


SUBROUTINE     INTEG3 

5    R  =  HI 
MMM  =  0 
K  =  KK 
AK  =  K 
IF     (JPRINT-2)     35.35.4400 

4400  PRINT  4401 

4401  FORMAT(50H  SUBROUTINE  INTE63  CALLED 
35  XP(K,KY)=(HI**K)*1.0E-4 

37  XPP(K,KY)=AK*(HI**(K-1) )*1.0E-4 
39  XM(K,KY) =XP(K,KY) 
41  XMP(K,KY)=XPP(K,KY) 

51  H=HI 

52  F=AMASS/(AMASS+1. 00898) 
54  FF  =  F 

IF  (MPRINT-1)  4090.4088.4090 

4088  PRINT  4089 

4089  FORMAT  ( 50H   WAVE  FUNCTION  VS  RADIUS 
PROINT=0 

IF(JPUNCH-l)  4090.4090,4087 
4087  PUNCH  4362 
4362  FORMAT (50H   ENERGY 

PUNCH  4361.  EN(KY) 
4361  F0RMAT1E13.5) 

4090  DO  4210  1  =  1  ,1000 
1  =  1 

CALL  KRR1 

4126  R=R+H 

4127  IF(NOMAX-l)  4130,4128,4130 

4128  IF  (R-0.2*RADIUS)  4130,4129,4129 

4129  IFIEN(KY)-V)  4220,4130,4130 

4130  IF(R-RSTOP+.001)  4140,4220,4220 
4140  IFIR-RI I+.001 )  4150,4160,4160 
4150  IF(R-RI+.001 )  4200,4180,4180 
4160  H=HIII 

4170  GO  TO  4210 

4180  H=HI I 

4190  GO  TO  4210 

4200  H=HI 

4210  CONTINUE 

4220  IF  (MSPIN)  4310,4230,4300 

43  00  PROINT=PROINT+(H/2.0)*XPLUS**2 

GO  TO  4320 
4310  PROINT=PROINT+ ( H/2 . 0 ) *XM I NUS**2 
4320  IF( MPRINT-1 142  30,4330, 4230 
4330  PRINT  4340 
4340  F0RMAT(5CH   NORMALIZATION 

4350  PRINT  436CPR0INT 
IF(JPUNCH-l)  4230,4230,4351 

4351  PUNCH  4340 
PUNCH  436CPROINT 

4360  F0RMAT(1E15.7) 
4230  RETURN 
END 
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Appendix  4.  Fortran  Statements  and  Sample  Output — continued 


5000  SUBROUTINE  KRR1 

ELL=K-1 

IF  (MSPIN)  5060,4125,5020 
5020  XPLUS=XP(K,KY) 

PROINT=PROINT+H*XPLUS**2 
5040  XPLUSP=XPP(K,KY) 

56  RR=R 

57  V=VPLUS( I ) 
59  X=XPLUS 

61  CALL  FR 

62  AI=FRR*H**2/2.0 

63  RR=R+H/2.0 

65  X  =  XPLUS+  XPLUSP*H/2.0+AI /4.0 
67  CALL  FR 

69  AI I  =FRR*(H*#2 1/2.0 

70  RR=R+H 

71  X=XPLUS+XPLUSP*H+AI I 
73  CALL  FR 

75  AI I I=FRR*H**2/2.0 
141  XPLUS=XPLUS+H*XPLUSP+(AI+2.0*AI I )/3.0 
143  XPLUSP  =XPLUSP  +(AI  +4 . 0*A I  I  +A I  I  I  ) / ( 3 . 0*H ) 
180  XP(K,KY)=XPLUS 
184  XPP( K,KY)=XPLUSP 
GO  TO  4110 
5060  XMINUS=XM(K,KY ) 

PROINT=PROINT+H*XMINUS**2 
5080  XMINUP=XMP(K,KY) 
77  RR=R 
98  X=XMINUS 

100  V=VMINUS( I ) 

101  CALL  FR 

103  CI=FRR*H**2/2.0 

105  RR=R+H/2.0 

107  X=XMINUS+XMINUP*H/2.0+CI/4.0 

109  CALL  FR 

111  CII=FRR*(H**2 ) /2.0 

112  RR=R+H 

113  X=XMINUS+XMINUP*H+CII 
115  CALL  FR 

117  CI  I  I  =FRR*H*#2/2.0 

149  XMlNUS  =  XMINUS  +  H*XMINUP+(CI+2.0-x-CI  I  )/3.0 
151  XMINUP=  XMINUP  + ( C I +4.0*CI I +C I  I  I  )  /  (  3 . 0*H ) 
188  XM(K,KY)=XMINUS 
192  XMP(K»KY)=XMINUP 
4110  MMM=MMM+1 

IF(JPRINT-l)  4114,4114,4120 
4120  I F( 1-3 )  4112,4112,4114 

4112  PRINT  4113,XP(K,KY)  ,XM(K,K,Y)  ,XPP(K,KY)  ,XMP(K,KY) 
PRINT  4113,XPLUS,XMINUS,XPLUSP,XMINUP 
PRINT  4113,AI,AII,AIII,CI,CII,CIII,H,R 
PRINT  4113,FRR,V,FF,X,RR 

4114  IF(MPRINT-l)  412  5,4115,4125 

4115  IF  (MSPIN)  4118,4118,4117 
4117  PS  I  (MMM)=XP(K,KY) 

RRR(MMM)=RR 
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Appendix  h.     Fortran  Statements  and  Sample  Output — continued 

GO  TO  4119 

4118  PSI (MMM)=XM(K,KY) 
RRR(MMM)=RR 

4119  IFIMMM-4)  4125,4122*412  5 
4113  FORMAT  (8E15.7) 

4122  PRINT  41 16* PS  I  (1)*RRR(1).PSI(2)»RRR(2)*PSI(3)»RRR(3)»PSI(4)  ,RRR{4) 
4116  FORMAT  ( E13.5 . F5 .2 . E13 . 5 ,F5 . 2 ,E1 3 . 5 , F5 . 2 , E13 . 5 . F5 . 2 ) 

IF( JPUNCH-1 )4135»4135,4132 
4132  PUNCH  4116»PSI(1)»RRR(1)»PSI(2)»RRR(2)»PSI(3)*RRR(3)»PSI(4)»RRR(4) 
4135  MMM=0 
4125  CONTINUE 

RETURN 

END 


SUBROUTINE  FR 
44  FRR=-.04826*FF*(EN(KY)-2  0.721*ELL*( ELL+1.0 )  / ( FF*RR**2 )-V ) *X 
•  RETURN 
END 
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Appendix  k.     Fortran  Statements  and  Sample  Output — continued 

SUBROUTINE    MATCH 
DO    6200    KY=1»KZ 
MPRINT=2 
KY  =  KY 

RSTOP=RSTART 
CALL     INTEG3 
PRINT    6180 
6180    F0RMATI119H  ENERGY  DIFFERENCE  FL     INSIDE  FL    OUT 

1SIDE       WAVEFNPLUS(R)       WAVEFNM I NUS ( R )  R  ) 

NTRY=1 
NO  =  l 

203  ENERGY(NO»KY)=EN(KY) 
IF    (MSPIN)    208.208.204 

204  AJLP(K,KY)=XP(K,KY)*XPP(K,KY> 
AILP(K.KY)=    XP(K.KY)**2 

REFL     (K,KY)=R*AJLP(K,KY)/AILP(K,KY) 

GO    TO    212 
2  08    AJLM(K,KY)=XM(K,KY)*XMP(K,KY) 
210    AILM(K.KY)=    XM(K»KY)**2 

REFL     (K,KY)=R*AJLM(K,KY)/AILM(K,KY) 
212    IF(MODE-l)     218.218.214 
214    DIF(NO)=REFL(KK.KY)-FLOUT(KY) 

218  IF(EN(KY))  219.219.220 

219  PMOM=SQRTF(-.04826*F*EN(KY)  ) 
IF  (MODE-1)  2191.2191.6001 

2191  XX(KY)=R*PMOM 

HTILDE(1,KY)=1.0+XX(KY) 

SPHTAN(1»KY)=-XX(KY) 

DO  6000  K=2»KK 

AK  =  K 

HTILDE(K,KY)=(2.0*AK-1.0)+(XX(KY)**2)/HTILDE(K-liKY) 
6000  SPHTAN(K,KY)=AK-HTILDE(K,KY) 

DIF(NO)=REFL(KK,KY)-SPHTAN(KK,KY) 

FLOUT (KY)=SPHTAN(KK,KY) 

GO  TO  6001 

220  PMOM  =     SQRTF( .04826*F*EN(KY) ) 
IF  (MODE-1)  2201.2201.6001 
XXIKY  )=R*PMOM 
S( l.KY)=XX(KY) 
DELTA(  1»KY)=0.0 
DO  226  K=1,KK 
AK  =  K 

S(K+l.KY)=(XX(KY)**2)*S(K,KY)/{ ( AK-DELTA ( K , KY ) ) **2+S ( K, KY ) **2 ) 
DELTA  (K+1»KY).  =  (AK-DELTA<K,KY)  )  *  (  S  (  K+l  ,KY  )  /S  (  K  ,  KY  )  )-AK 
P(1.KY)=1.0 
P(2.KY)=1.0/XX(KY) 
Q( l.KY)=0.0 
Q(2»KY)=-1.0 
DO  254  K=3»KK 
AK  =  K 

P(K.KY)=(2.0*AK-3.0)*P(K-1,KY)/XX(KY)-P( K-2.KY) 
Q(K,KY )=(2.0*AK-3.0)*Q(K-1 .KY ) /XX ( KY ) -Q( K-2 »KY ) 
DO  268  K=1.KK 
GMF (K,KY)=(P(K,KY) **2~Q ( K , KY ) **2 ) *COSF ( 2 . 0*XX ( KY > )~2 .0*P ( K ,KY ) * 

1Q(K,KY)*SINF(2.0*XX(KY) ) 


2201 

221 

222 

223 

224 

226 

240 

244 

246 

248 

250 

251 

252 

254 

255 

256 
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Appendix  h.     Fortran  Statements  and  Sample  Output--uontinued 

2  58  TFG(K,KY)=2.0*P(K,KY)*Q(K,KY)*COSF(2.0*XX(KY)  )  +  ( P ( K,KY ) **2 

1-Q(K,KY)**2)*SINF(2.0*XX(KY)  ) 
268  FPG(K,KY)=P(K,KY)**2+Q(K»KY)**2 

FLOUT(KY)=DELTA(KK,KY)-S(KK,KY)*TFG(KK,KY)/(GMF(KK,KY)+FPG(KK,KY) ) 

DIF(NO)=REFL(KK,KY)-FLOUT(KY) 

IF(JPRINT-2)  6001,6001,271 

271  PRINT  272,XX(KY) ,S(KK,KY) ,DELTA(KK,KY) ,P(KK,KY) ,Q(KK,KY) , 
1GMF(KK,KY) ,TFG(KK,KY) ,FPG(KK,KY) 

272  FORMAK8E15.7) 

6001  IF(NTRY-2)    6002,6002,6075 

6002  IF(NO-NOMAX)    6004,6100,6100 

6004  IF(NO-2)    6005,6080,6090 

6005  IF(DIF(NO)+1.0)  6010,6010,6020 
6010  IF(DIFINO)-l.O) '  6030,6040,6040 
6020    EN(KY)=( 1.0-CONV)*EN(KY) 

GO    TO    6100 
6030    EN(KY  )  =  ( 1.0+CONV )*EN ( KY ) 

GO    TO    6100 
6040    IF(DIF(NO)+0.1  )  6050,6050,6060 

6050    IF(DIF(NO)-0.1)  6070,6100,6100 

6060    EN(KY)=( 1 .0-0.2 *CONV ) *EN ( KY ) 

GO    TO    6100 
60  70    EN(KY)=(1.0+0.2*CONV)*EN(KY) 

GO    TO    6100 

6075  IFJNO-2)  6076,6080,6090 

6076  EN(KY)  =  (E(NTRY-l)-E(NTRY-2)  )  /S IGMAE+E ( NTRY-1 ) 
GO  TO  6100 

608  0  EN (KY )=( ENERGY ( NO-1 »KY)*DIF (NO) -ENERGY ( NO »KY)*D I F(NO-l) ) / 

1(DIF(N0)-DIF(N0-1) ) 
GO  TO  6100 
6090  EN(KY)=( ( ( ENERGY ( N0-2 »KY ) *D I F ( NO-1 ) -ENERGY ( NO- 1 ,KY ) *D I F ( N0-2 ) )* 

lDIF(N0)/(DIF(N0-l)-DIF(N0-2 ) ) )-( ( ENERGY ( NO-1 ,KY ) *DI F ( NO ) 

2-ENERGY(N0,KY)*DIF(N0-l )  )*DIF(N0-2 ) / ( DIF ( NO )-DI F ( NO-l )  )  )  )/ 

3(DIF(N0)-DIF(N0-2) ) 
6100  PRINT  6160, ENERGY! NO, KY)  ,DIF(N0) ,REFL(KK,KY) , FLOUT (KY)  ,XP(KK,KY)  , 

1XM(KK,KY)  ,R 
6160  FORMAT(F15.7,5E15.7»F8.3) 

IF  (NO-NOMAX)  6114,6200,6200 

6114  IFIDIF  (N0)+  EPS)  6120,6120,6115 

6115  IFIDIF  (NO)-  EPS)  6130,6120,6120 

6120  IF  (JPRINT-1)  6123,6123,6121 

6121  PRINT  6160,  EN(KY) 
6123  CALL  INTEG3 

NO=NO+l 
GO  TO  203 

6130  E(NTRY)=ENERGY(NO,KY) 
IFfRSTOP  -RMAX)  6131,6132,6132 

6131  IF(JPRINT-2)6150, 6132, 6132 

6132  MPRINT=1 
6140  CALL  INTEG3 
6150  MPRINT=2 

IF  (RSTOP-RMAX)  6165,6200,6200 
6165  RSTOP=RSTOP+LOGF(SIGMAR) /PMOM 
IF  (RSTOP-RMAX)  6168,6168,6167 

6167  RSTOP=RMAX 

6168  NTRY=NTRY+1 

CALL  INTEG3 
NO=l 

GO  TO  203 
6200  CONTINUE 
RETURN 
END 
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